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SMALL  SCALE  STRUCTURE  IN  THE  EARTH’S  IONOSPHERE: 

THEORY  AND  NUMERICAL  SIMULATION 

1.  Introduction 

It  Is  generally  believed  that  the  existence  of  Ionospheric  structures 
with  scale  sizes  of  tens  of  kilometers  or  smaller  can  be  attributed 
primarily  to  the  onset  and  evolution  of  instabilities  of  one  sort  or 
another.  These  Instabilities  can  be  thought  of  as  being  superimposed  on, 
and  indeed  evolving  from,  Che  larger  scale  ionospheric  configuration. 
Among  the  numerous  such  structures  it  is  usually  only  those  that  are  of 
reasonably  large  amplitude  or  those  which  cause  problems  (e.g. 
communications  interference)  that  attract  interest  and  study.  Still,  this 
number  is  greater  than  we  can  reasonably  treat  here.  We  shall  therefore 
limit  our  discussion  to  two  such  structures  whose  physics  and  evolution  we 
believe  we  understand  reasonably  well:  1)  the  steepening  and  subsequent 
recursive  splitting  of  barium  clouds  released  in  the  ionosphere,  driven  by 
the  gradient  drift  instability;  and  2)  the  formation  and  buoyant  rise  of 
low  density  "bubbles"  of  plasma  in  the  nighttime  equatorial  ionosphere, 
known  as  equatorial  spread  7  (ESF),  driven  by  the  colllslonal  Rayleigh- 
Taylor  instability.  Each  of  these  instabilities  derive  from  the  same  set 
of  plasma  fluid  equations  and  the  same  set  of  physical  approximations, 
differing  only  in  geometry  and  in  the  identity  of  the  driving  terms;  hence 
we  shall  attempt  here  to  unify  their  description  as  much  as  possible.  We 
shall  find  that  one  of  the  characteristics  of  structures  resulting  from 
chese  instabilities  is  their  tendency  to  be  "field  aligned",  that  is,  for 
the  plasma  gradients  and  velocities  parallel  to  the  magnetic  field  to  be 
much  smaller  than  those  perpendicular  to  the  magnetic  field.  Our 
discussion  will  therefore  focus  on  plasma  motion  perpendicular  to  the 
ambient  magnetic  field. 

In  Figure  1  we  show  a  photograph  of  the  Spruce  event,  a  barium  cloud 
released  at  188  km  altitude  in  February  of  1971,  24  minutes  after 

release.  The  cloud  was  originally  released  as  a  gaussian  distribution  of 
neutral  barium  which  was  subsequently  photolonized  by  sunlight.  In  the 
very  center  of  the  photograph,  our  line  of  sight  is  parallel  to  the 
magnetic  field  lines  at  the  cloud  altitude,  revealing  the  fine  scale 
structure  (termed  "strlations")  that  has  evolved  from  this  originally 
Manuacript  approved  April  27,  1983. 
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nearly  gaussian  distribution  of  plasma.  In  Figure  2  we  show  a  sketch  of 
what  we  believe  to  be  the  typical  evolution  of  a  barium  cloud  like  Spruce, 
dervived  from  experimental  observations  and  numerical  simulations.  The 
inital  steepening  of  the  top  of  the  two-dimensinal  cloud  is  caused  by  the 
buildup  of  polarization  charge  on  its  sides,  causing  the  high  density 
center  of  the  cloud  to  |  x  B  drift  in  the  direction  of  the  neutral  wind  to 
a  greater  extent  than  the  low  density  edges.  As  the  plasma  gradient  on  the 
top  of  the  cloud  becomes  steeper,  the  growth  rate  of  the  gradient  drift 
instability  (to  be  described  later)  active  there  becomes  larger  and 
eventually  small  perturbations  on  this  gradient  are  amplified  into  visible 
ripples,  which  in  turn  evolve  into  finger-like  structures.  Each  of  the 
strucutes  emerging  from  the  steepened  edge  of  the  cloud  then  evolve  into 
smaller  clouds,  and  the  process  begins  again,  resulting  in  a  cascade  of 
recursively  decreasing  scale  sizes  until  the  instability  is  stopped  by 
dissipation  or  other  mechanisms  which  act  more  effectively  on  the  smaller 
space  scales. 

In  Figure  3  we  show  maps  of  1-m  irregularities  taken  from  Tsunoda 
(1981)  at  the  earth's  magnetic  equator  during  equatorial  spread  F  (ESF). 
These  Irregularities  have  been  shown  to  be  closely  associated  with 
"bubbles"  or  regions  of  large  electron  density  depletion  in  the  equatorial 
ionosphere,  and  can  be  thought  of  as  at  least  a  partial  map  of  the 
locations  of  severe  electron  density  depletion.  In  Figure  4  we  show  the 
results  of  a  numerical  simulation  from  Zalesak,  et  al.  (1982),  showing  the 
time  evolution  of  electron  density  contours  at  the  earth's  equator.  The 
equatorial  ionosophere  was  originally  laminar  with  a  maximum  in  electron 
density  at  430  km  altitude.  A  sinusoidal  perturbation  was  applied  in  the 
east-west  direction.  The  results  show  that  the  observed  "bubbles"  consist 
of  low  density  plasma  which  has  been  transported  from  very  low  altitudes  up 
through  the  F2  peak  and  beyond  by  the  nonlinear  evolution  of  the 
colllsional  Raylelgh-Taylor  instability.  The  westward  and  eastward  tilts 
of  the  bubble  are  due  to  an  eastward  neutral  wind  blowing  at  the  equator 
coupled  along  magnetic  field  lines  to  background  ionization  (e.g.,  E 
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regions)  at  higher  and  lover  latitudes.  Note  that  the  various  tilts  of  the 
bubbles  In  Figures  3  and  4  are  consistent  when  allowance  Is  made  for  the 
reversed  abscissae  In  the  two  plots. 

In  Section  2  we  shall  present  a  qualitive,  physical  description  of  the 
instabilities  active  in  ionospheric  barium  cloud  and  equatorial  spread  F 
(ESF)  cases.  In  Section  3  we  derive  the  set  of  equations  describing  the 
motion  of  ionospheric  plasma  in  general,  and  the  evolution  of  barium  clouds 
and  ESF  bubbles  in  particular.  In  Section  4  we  discuss  the  simplifications 
made  in  constructing  a  mathematical  representation  of  the  physical 
system.  In  Section  5  we  derive  and  summarize  the  equations  describing  the 
"simplest  case"  geometries  and  assumptions  for  each  of  the  instabilities. 
Finally,  in  Sections  6  through  8,  we  treat  the  numerical  integration  of 
these  differential  equations. 

2.  The  Gradient  Drif t/Colllsional  Rayleigh-Taylor  Instability 

In  this  section  we  shall  attempt  to  give  a  qualitative  physical 
picture  of  the  gradient  drift  and  collional  Rayleigh-Taylor  Instabilities, 
both  of  which  are  caused  by  the  differential  motion  of  ions  and  electrons 
perpendicular  to  the  magnetic  field.  We  consider  a  two-dimensional  x-y 
plane  perpendicular  to  the  ambient  magnetic  field  B.  A  plasma  species  a  in 
this  plane  embedded  in  a  neutral  gas  will  respond  to  an  external  force 
perpendicular  to  g,  in  two  ways:  1)  by  drifting  in  a  direction 
perpendicular  to  both  g  and  F,^  (Hall  mobility)  and  2)  usually  to  a 
lesser  extent,  by  drifting  in  a  direction  parallel  to  F^  (Pedersen 
mobility).  We  shall  explicity  derive  these  drifts  in  Section  3.  We  shall 
take  our  plasma  to  consist  of  a  single  ion  species,  denoted  by  subscript  i, 
and  of  electrons,  denoted  by  subscript  e.  The  instabilities  under 
discussion  result  from  the  fact  that  the  ions  and  electrons  drift  with 
different  velocities  and  directions  in  response  to  the  same  external 
force.  In  regions  where  plasma  density  gradients  exist,  this  difference  in 
velocities  causes  polarization  charge  to  be  created  in  our  originally 


neutral  plasma,  which  in  turn  produce*  a  polarisation  electric  field.  The 
plaeaa  drift  associated  with  the  electric  field  will  cause  growth  of  a 
perturbation  when  the  plasaa  gradient  is  properly  aligned. 

In  Figure  5  we  show  contours  of  content  plasaa  density  n  in  the  two- 
dlaenslonal  x-y  plane,  where  we  assuae  the  aagnetlc  field  g,  to  be  aligned 
along  the  positive  s  axis.  Depicted  is  a  one-dimensional  "slab**  of  plasaa 
n(y)  such  that  n  aaxlalzes  at  y*yQ,  super laposed  on  which  is  a  sinusoidal 
perturbation  proportional  to  sin  kx,  where  k  is  a  vavenuaber.  Either  a 
downward-directed  gravational  acceleration  (in  the  collislonal  Rayleigh- 
Taylor  instability)  or  a  downward-directed  neutral  wind  (in  the  gradient 
drift  instability)  will  cause  the  ions  to  drift  leftward  relative  to  the 
electrons,  leaving  polarization  charge  where  the  relative  drift  has 
components  parallel  to  the  density  gradient,  as  indicated  in  Figure  5. 
This  polarization  charge  induces  a  polarization  electric  field  which  in 
turn  induces  an  additional  plasma  drift  in  the  Ep  x  B  direction.  This 
drift  is  such  as  to  enhance  the  perturbation  for  y  <  yQ  (instability),  as 
seen  in  Figure  5.  In  their  most  simplified  geometries  the  linear  growth 
rates  y  for  the  gradient  drift  and  collislonal  Rayleigh-Taylor 
instabilities  are 

y  -  -  ~  (gradient  drift)  (1) 

Y  ■  “  g  *  (Rayleigh-Taylor)  (2) 

where  u^n  is  the  ion-neutral  collision  frequency.  We  note  here  that  the 
gradient  drift  instability  may  be  thought  of  as  being  driven  by  an  ambient 
electric  field  simply  by  performing  a  Lorentz  transformation  into  the  rest 
frame  of  the  neutral  gas. 

The  above  picture  of  the  early  (linear)  stage  of  the  instability 
evolution  is  quite  informative,  but  unfortunately  falls  short  of 
illuminating  the  complex  nonlinear  evolution  of  barium  clouds  and 
equatorial  spread  F  "bubbles”.  In  the  next  section  we  derive  the  equations 


necessary  for  s  complete  nonlinear  description  of  these  phenomena,  vhlch  In 
general  require  numerical  techniques  for  their  solution. 

3.  The  Motion  of  Ionospheric  Plasma 

We  shall  be  concerned  here  with  the  motion  of  plasma  consisting  of  ions  and 
electrons  In  the  presence  of  a  neutral  gas  and  magnetic  field  B,  subject  to 
an  external  force.  We  shall  also  be  interested  in  the  electric 
current  J  arising  from  the  differential  motion  of  thevarlous  species 
comprising  the  plasma.  In  the  course  of  deriving  the  equations  we  shall 
make  some  assumptions  which  are  crucial  to  the  model: 

1)  We  shall  assume  the  plasma  can  be  adequately  described  by  the 
fluid  approximation.  This  assumes  that  the  effective  collision  rate  of 
each  plasma  species  with  itself  is  sufficiently  high  to  maintain  near 
Maxwellian  distribution  functions  on  time  scales  short  compared  to  the 
times  of  interest,  and  is  well  satisfied  for  the  plasmas  we  treat  here. 

2)  We  shall  assume  that  the  electric  fields  E  are  electrostatic 
(i.e.,  V  x  E  *  0)  and  hence  can  be  described  using  a  scalar  potential  <j> 
such  that  E  »  -  7  $.  Note  that  this  implies  3B/at  *  0.  The  validity  of 
this  assumption  can  be  related  to  the  fact  that  the  Alfven  velocity  is  much 
larger  than  any  other  propagation  speeds  of  interest  for  the  plasmas  we 
treat  here.  The  assumption  is  also  checked  a  posteriori  by  verifying  that 
the  calculated  currents  and  displacement  currents  produce  negligible  time 
variations  in  ft  which  in  turn  produce  negligible  7  x  E. 

3)  We  assume  plasma  quasi-neutrality;  that  is, 

l  nt  qi  -  Vs  (3) 

where  n  is  number  density,  q  is  ion  specie'  charge,  e  is  the  electron 
charge,  the  subscripts  i  and  e  refer  to  ions  and  electrons  respectively, 
and  the  sum  is  taken  over  all  ion  species.  This  assumption  is  a  statement 
that  the  Debye  length  is  small  compared  to  all  length  scales  of  interest, 
and  again  can  be  verified  a  posteriori  by  evaluating  7  •  E.  Note  that 
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this  assumption  laplies  that  7  •  J  ■  0,  where  J  is  the  electric  current. 
In  addition  to  the  above  there  are  some  other  aaauaptions  which,  while  they 
are  not  essential  to  the  basic  model,  are  nonetheless  valid  for  many  of  the 
physical  situations  which  we  shall  treat  and  Impart  a  simplicity  which  we 
shall  find  convenient  here: 

4)  We  shall  assume  the  electrostatic  potential  $  to  be  constant  along 
magnetic  field  lines.  As  we  shall  see  later,  the  electrical  conductivity 
along  magnetic  field  lines  is  much  greater  that  that  perpendicular  to 
magnetic  field  lines,  meaning  that  appreciable  differences  in  potential 
along  a  field  line  will  quickly  be  reduced  by  the  resultant  current.  This 
assumption  will  break  down  for  sufficiently  small  scale  lengths 
perpendicular  to  the  magnetic  field,  and  for  sufficiently  large  distances 
along  the  magnetic  field. 

5)  We  shall  assume  that  the  inertial  terms  in  the  plasma  species 
momentum  equations,  i.e.,  the  left  hand  side  of  Equation  <5),  are 
negligible  with  respect  to  the  other  terms  in  the  equation.  This 
assumption  is  justified  whenever  the  time  scales  of  interest  are  longer 
than  the  mean  time  between  collisions  for  ions. 

6)  We  shall  neglect  all  collisions  between  species  except  those 
between  ions  and  the  neutral  gas.  This  is  justified  simply  by  an  evalua¬ 
tion  of  the  magnitudes  of  the  terms  involved. 

7)  We  shall  ignore  production  and  loss  terms  which  may  appear  as 
sources  and  sinks  in  the  plasma  continuity  equations  as  a  result  of 
chemistry,  photoionization,  etc. 

Assumptions  (4)  through  (7)  above,  although  made  in  this  paper,  are 
not  necessary  within  the  theoretical  and  computational  framework  we  have 
developed,  and  adequate  means  exist  to  delete  them,  if  necessary. 

The  continuity  and  momentum  equations  describing  plasma  species  a  are: 
v  .  (n  vp  -  0  (4) 


6 


q_  v  *  B 

7)  v  -  -S  (E  +  ^L_Z) 

J  cl  n  c  * 


vln  <1.  ’  V  -*■■ 


n  ■ 

a  a 


i 


(5) 


where  Che  subscript  a  denotes  the  plasma  species  (i  for  ions,  e  for 
electrons,  for  example),  n  is  the  species  number  density,  v  is  the  species 
fluid  velocity,  P  is  pressure,  g,  is  the  electric  field,  g  is  the 
gravitational  acceleration,  q  is  the  species  charge,  is  the  species 
collision  frequency  with  the  neutral  gas,  U  is  the  neutral  wind  velocity, 
c  is  the  speed  of  light,  and  m  is  the  species  particle  mass.  He  can 
rewrite  this  equation  as 


F  /m 
~a  a 


v  v  • 
an  ~a 


0 


where 


(6) 


F  =  q  E  +  mg  +  v„.m  U  -  7P  /n 
~a  a  ~  a  an  a  m i  n  a 


~  (4r  +  v  •  7)  v  m 
a  t  ~a  ~a  a 


(7) 


If  we  place  ourselves  in  a  Cartesian  coordinate  system  in  which  £  is 
aligned  along  the  z  axis,  and  if  we  treat  as  a  given  quantity  then  a 
componentwise  evaluation  of  Equation  (6)  yields  a  set  of  three  equations  in 
three  unknowns,  the  three  components  of  v  .  The  formal  solution  is 


v  i" 

~olL 


k.  F  +  k,  F  x 
lot  ~al  2a  ~al 


(8) 


~al 


k  F. 
oa  ~l 


(9) 


where 
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k  -  VqP  c  fi  _  i 

la  naF7T[  l  +  ^an,fl.,"] 

(u  /n  )2 

k2a  'ql  t1  ‘  1  ^ (Va  7Q  )2^ 

a  an  a 


oa 


a  an 


-1 


(10) 

(11) 

(12) 


z  s  ft/lBl 


(13) 


(14) 


The  vector  subscripts  1  and  I  refer  to  the  components  of  the  vector  which 
are  perpendicular  and  parallel  respectively  to  z.  The 
quantities  k^,  k£,  and  k0  above  are  referred  to  as  the  Pedersen,  Hall,  and 
direct  mobilities  respectively.  It  should  be  pointed  out  that  Equations 
(8)  and  (9)  are  only  truly  closed  form  expressions  when  the  Inertial  terms 
(the  last  term  on  the  right  hand  side  of  Equation  (7))  are  neglected,  an 
assumption  we  have  made  previously.  Typical  ranges  for  collision 
frequencies  are:  ~  30  sec"*,  ven  -  800  sec"*  at  150  km  altitude; 
and  v^n  ~  10“*sec"*,  ven  ~  1  sec"*  at  500  km  altitude. 

As  we  will  see  later,  we  will  use  the  concept  of  "layers”  to 
distinguish  the  various  ion  species,  so  for  the  moment  we  can  consider  only 
a  single  ion  species,  denoted  by  subscript  i,  and  the  associated  electrons, 
denoted  by  subscript  e.  We  will  also  consider  only  singly  charged  ions  so 
that  q^  •  e  and  qg  ■  -  e.  Noting  that  ven/^e  *  0  we  obtain 


kli 

.  in  R  c 

—  ri7TbT 

(15) 

kle 

-  0 

(16) 

k2i 

m  hiw 

(17) 

k2e 

"  " 

(18) 
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where 


Rt  2  (1  +  Vin2/^2)-1  (19) 

We  now  define  Che  perpendicular  current 

J  =  V  n  q  v  (20) 

~1  L  a  a  ad 
a 

Substituting  Equations  (15)  through  (18)  and  (8)  Into  Equation  (20),  and 
using  the  quasi-neutrality  approximation 


nl 


n  =  n 
e 


we  obtain 


(21) 


in  nc 

~i  n1  i  TbT  ~ii 


+  T-  (Ri  hi  +  lel)x  2 


(22) 


For  the  barium  cloud  and  equatorial  spread  F  (ESF)  problems  we  shall  treat 
here,  we  shall  only  consider  neutral  winds,  electric  fields,  and  gravity  as 
external  forces.  Hence 


F . .  ■  e  E .  +  m,  g ,  +  v.  m.  U  . 

~il  ~1  i  ftl  in  i  ~nl 

F  *  -  e  E  .  +  m  g. 

~ei  ~1  e  *-1 

Note  that  we  have  neglected  the  small  term  ven  mg  in  Equation  (24). 
We  obtain 


(23) 

(24) 
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i  fir  <e  ~1  +  "i  ftl  +  vm  “i  Hni> 


J 


1 


R 


(25) 


+  R. 


if  +  Si 


+  v 


in 


m  U 
i  ~nl 


X  Z 


Since  0.01  <  <  1.0  we  may  neglect  nig/Rj  with  respect  to  n^. 

Defining  the  Pedersen  conductivity 

v. 

in  nee 

°p  =  Ri  ni  TIT 

and  noting  that  1  -  we  obtain 


(26) 


J 


1 


m  m 

a  [e  +  —  %  +  v,  "■  U 

p  e  *1  in  e  ~ni 


v.  n.  a.  m. 

+  (  -  -il  E  +  -i-i  g  +  Q.  — 
2.  ~1  v.  e  i  e 

i  in 


U  lx  zl 
~nl'  J 


(27) 


Our  need  for  an  expression  for  stems  from  our  need  for  its  divergence  to 
evaluate  7  •  £  (  »  0  by  quasi-neutrality),  as  we  shall  see  in  the  next 
section. 


4.  Model  Simplification  and  Mathematical  Representation 

We  shall  model  our  physical  system  using  a  simplified  model  as 
depicted  in  Figure  6.  The  magnetic  field  lines  are  assumed  to  be  straight, 
to  be  aligned  along  the  z  axis  of  our  cartesian  coordinate  system,  and  to 
terminate  in  Insulators  at  z  ■  +_  °».  The  plasma  of  interest  is  threaded  by 
these  magnetic  field  lines,  and  is  divided  into  thin  planes  or  ’’layers'’  of 
plasma  perpendicular  to  the  magnetic  field.  Since  we  have  neglected 
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collisions  between  different  plssas  species,  we  nay  use  the  device  of 
layers  to  treat  multiple  ion  species  at  a  single  point  in  space  sinply  by 
allowing  multiple  layers  to  occupy  the  same  plane  in  space,  one  for  each 
ion  species.  In  this  way  a  “layer”  consists  only  of  a  single  ion  species 
and  its  associated  electrons. 

Our  quasi-neutrality  assumption  demands  that 


V 


J 


1“  J 
3x  x 


J 

y 


j 

z 


0 


(28) 


Integrating  Equation  (28)  along  z  and  noting  from  Figure  6  that  Jz  vanishes 
at  z  *  +  •  we  obtain 


/  +  "  V  *  Ji  dz  *  0 


(29) 


where 


V 


1 


(30) 


From  our  model  as  depicted  in  Figure  6  we  may  approximate  the  Integral  in 
Equation  (29)  by  a  discrete  sum 


7i  •  ilk  “k  ■  0  <31) 

where  the  subscript  k  refers  to  the  layer  number,  N  is  the  total  number  of 
layers  in  the  system,  and  Azk  is  the  thickness  of  layer  k  measured  along 
the  magnetic  field  line.  By  our  assumption  of  equlpotential  magnetic  field 
lines  and  electrostatic  electric  fields 

Eik(*»y)  ■  "  7x<!>(x,y)  for  all  k  (32) 
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aaatfMiniMi 


Then  Equation  (31)  becomes 


(33) 


!_  r^ia  j  ii)  +L  r^a  e  ii) 

3x  '•fl  p  3yj  k  3y  '■Q  p  3x;  k 


M  i_  r^4a  s  )  +  n  l_  rlia  s  ) 

3y  3x  (jj  V  k  3x  3y  ^  pJ  k 


Text  _  r°i  mi  . 

~lk  "  ^pk  ( e  ®1  Vin  e  Mil 

»N/ 

n.m  m 

+  (“=“^  8,  +  fl.  —  Ujx  *]. 

'v  e  1  1  e  Jk 

in  ~ 


and  the  subscript  k  denoting  layer  number  on  terms  within  parenthesis 
operates  on  all  terms  within  those  parentheses.  Equation  (33)  is  a  second 
order  elliptic  partial  differential  equation  for  $(x,y),  subject  to 
boundary  conditions  on  Our  reason  for  writing  Equation  (33)  in  the  form 
we  did  is  related  to  the  following  picture  of  the  physics.  The  external 
forces  acting  on  a  plasma,  in  this  case  gravity  and  a  neutral  wind 
collision  term,  will  Induce  a  current  to  flow.  In  general  this  current 
will  not  satisfy  7  »  J  •  0,  meaning  that  in  certain  regions  there  will  be  a 
build-up  of  polarization  charge,  resulting  in  an  electric  field  which 
causes  secondary  currents  to  flow.  Over  time  scales  much  shorter  than 
those  of  Interest  here,  a  quasi-steady  state  is  reached  such  that 
subsequent  plasma  motion  is  well  described  by  7  •  ^  ■  0.  In  this  physical 
picture  the  electric  field  represents  the  response  of  the  plasma  to  a  given 
externally  driven  current  system.  Thus  the  right  hand  side  of  the  Equation 
(33)  may  be  regarded  as  the  "known"  divergence  of  the  external  current, 
which  we  shall  denote  below  by  R,  and  the  left  hand  side  regarded  as  a 
differential  operator  L  operating  on 
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-  R  (37) 

The  operator  L  li  a  hermit lan  operator  in  the  limit  as  the  "Hall  terms" 
may  be  neglected,  as  is  often  the  case  at  higher  altitudes  where  is 

small. 


5.  The  simplest  Case  Equations  for  Barium  Clouds  and  for  ESF 

The  simplest  case  equations  for  each  of  our  physical  systems  are  for  one 
level  only,  i.e.,  N-l,  and  for  altitudes  such  that  terms  of 
order  (\>in/Qi)2  may  be  neglected  with  respect  to  terms  of  order  (vin/ni). 
Ue  also  treat  only  one  external  force  for  each  case,  and  align  that  force 
along  one  of  the  coordinate  axes.  Since  we  have  only  one  level,  we  drop 
the  subscript  k. 

For  barium  clouds,  we  assume  B  to  be  aligned  along  the  z  axis,  that 

the  only  external  force  is  a  neutral  wind  U  =  U  7. 

~n  n 

Then 


Z 

P 


y  +  Q, 


A 

X) 


(38) 


Since  £p  is  already  of  order  (v^/fl^),  we  nay  neglect  the  first  term  with 
respect  to  the  second.  Then 


_  ,  BU 

,,  .  J*«  .»  a 
1  ~1  3x  p  c 


(39) 


where  we  have  used  Equation  (14). 

Noting  that  H  in  Equation  (33)  is  of  order  (vin/^i)2  we  obtain 


BU 

(Z  7.$)  -  (i  — -) 
p  3x  v  p  c  ' 


(40) 


For  the  equatorial  spread  F  case  we  assume  a  single  plane  of  plasma 

A 

located  at  the  magnetic  equator  such  that  B  is  along  the  z  axis  and  y  is 
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r 


"up".  Our  only  external  fore*  is  gravity  •  -  g?(g  ■  +  980  ca/sec2). 
Then 


(41) 


The  first  tera  is  of  order  vln/fl^  tiaes  the  second  and  aay  therefore  be 
neglected.  Then 


V 


1 


Again  neglecting  H  in  Equation  (33)  we  obtain 


(42) 


V 


1 


«p  V> 


(43) 


For  both  the  one-layer  bariua  cloud  and  ESF  cases,  one  aay  solve 
either  the  electron  or  the  ion  continuity  equations,  since  quasi-neutrality 
makes  them  equivalent  (but  not  Identical).  For  simplicity  we  choose  the 
electron  equation  since  we  may  neglect  the  Pedersen  terms 
there  (ven/Qe  -  0). 

Summarizing  the  equations  we  must  solve  for  each  case  we  get 


<ne  X.i>  ’ 


(44) 


71  #  UP  *  3S/3X 


(45) 


^el 


-  -  ~  F  x  z 
eB  ~el 


(46) 


F  . 
-el 


eV^$  for  bariua  clouds 

A 

e  -  m^gy  for  ESF 


(47) 
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s 


Ep  ^t/c  for  barium  clouds 
-  Ip  Bg/(cvlQ)  for  ESF 


(48) 


£  -  As  (v.  /Q. )nce/B  (49) 

p  in  1 

Solution  of  these  equations  requires  the  use  of  two-dimensional  numerical 
3 inula t ion  techniques. 


6.  Numerical  Simulation:  General 

We  saw  In  the  previous  section  that  In  the  simplest  case  for  the 
barium  cloud  and  equatorial  spread  F  (ESF)  problems,  we  can  reduce  our 
system  to  two  partial  differential  equations  posed  on  a  two  dimensional 
plane : 


It  +  \  •  <"  i.i>  ■ 0  <50> 

71  *  (Ep  71  **  "  3S/a*  (51) 

where  £p  and  S  are  explicity  given  functions  of  n  and  v^  is  an  explicity 
given  function  of  Equation  (50)  is  hyperbolic  while  Equation  (51)  Is 
elliptic.  Both  require  the  imposition  of  physically  relevant  boundary 
conditions.  Conceptually  one  solves  this  coupled  system  of  equations  as 
follows.  At  any  given  time  t,  we  assume  that  we  know  n(x,y,t)  and 
therefore  £p(x,y,t)  and  S(x,y,t).  We  can  then  solve  Equation  (51)  for  Its 
single  scalar  unknown  $(x,y,t),  given  properly  specified  boundary 
conditions  on  $  and/or  its  derivatives.  Knowing  $  we  can 

compute  ^^(x.y.t)  explicity.  We  can  then  solve  Equation  (50) 

for  n(x,y,t  +  At)  where  AC  is  a  small  time  increment.  The  process  Is 
repeated  recursively  until  the  solution  is  advanced  to  the  desired  time. 
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Within  the  above  context  several  numerical  approaches  exist  for 
solving  this  system  of  coupled  partial  differential  equations:  spectral 
methods,  finite  element  methods,  Galerkln  methods,  and  finite  difference 
methods,  among  others.  We  have  chosen  finite  difference  methods  here  for 
reasons  of  simplicity,  computational  efficiency,  and  most  Importantly 
because  acceptable  techniques  for  solving  Equation  (SO)  in  the  presence  of 
large  gradients  in  n  presently  exist  only  within  the  finite  difference 
domain.  Fundamental  to  finite  difference  techniques  is  their  use  of  a 
"grid",  that  is,  a  discrete  set  of  points  in  space  and  time  denoted 
by  (x^,y^,t°),  1  <  i  <  NX,  1  <  j  <  NY,  1  <  n  <  <■  where  i,  j,  m,  NX  and  NY 
are  integers,  on  which  the  solution  is  computed.  For  instance,  the 
electrostatic  potential  $(x,y,t)  at  x  •  x^,  y  ■  y^,  and 

t  ■  t®  would  be  denoted  In  Figure  7,  we  show  an  example  of  a  finite 
difference  grid  in  space,  and  we  also  show  how  the  grid  would  look  in  the 
case  of  multiple  layers  of  plasma,  although  we  shall  treat  only  a  single 
layer  here.  Note  that  the  four  "nearest  neighbors"  of  the  grid 
point  (xif  yj)  are  the  grid  points  (xi+i,  yj),  (xj-i,  yj), 

(Xi*y j+l>*  and  (xi»  yj**l)‘  Many  finite  difference  techniques  employ  what 
is  known  as  a  staggered  grid,  meaning  that  different  dependent  variables  (n 
and  <J>  for  instance)  are  evaluated  on  different  grids  in  space  and  possibly 
time.  We  do  not  employ  staggered  grids  here;  all  dependent  variables  are 
evaluated  on  exactly  the  same  grid. 

Looking  at  Equations  (50)  and  (51)  we  see  that  there  is  more  to 
distinguish  them  than  just  their  hyperbolicity  and  ellipticity 
respectively.  Both  equations  require  the  evaluation  of  spatial 
derivatives;  but  Equation  (50)  requires  in  addition  the  evaluation  of 
temporal  derivatives.  Precisely  because  we  do  not  yet  know  the  solution  at 
a  time  later  than  it  has  thus  far  been  computed,  the  treatment  of  temporal 
derivatives  is  qualitatively  different  from  that  of  spatial  derivatives. 
More  importantly,  it  has  been  found  empirically  that  it  is  the  numerical 
treatment  of  Equation  (50)  which  will  "make  or  break”  the  solution  to  the 
total  system  of  equations.  Specifically,  Equation  (51),  once  properly 
discretized  (i.e.,  once  the  spatial  derivatives  are  properly  represented  in 
finite  difference  form)  simply  yields  a  system  of  linear  equations,  albeit 
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a  very  large  system.  Our  experience  has  been  that  a  number  of  algorithms 
will  successfully  yield  a  solution  to  this  linear  system,  although  it  may 
be  difficult  to  find  one  algorithm  that  will  sclve  the  system  for  all 
possible  physical  parameters.  Accordingly  we  shall  discuss  the  numerical 
treatment  of  Equation  (51)  only  briefly  here,  in  the  next  section,  and 
reserve  the  bulk  of  our  discussion  for  Equation  (50). 


7.  The  Numerical  Solution  of  the  Potential  Equation 

As  was  mentioned  in  the  previous  section,  the  numerical  solution  of 
Equation  (51)  takes  place  in  two  stages:  1)  the  discretization  of  the 

spatial  derivatives  and  boundary  conditions  in  finite  difference  form, 
resulting  in  a  large  linear  system  of  NX  •  NY  equations  for  the 
NX  •  NY  unknowns  and  2)  the  solution  of  this  large  linear  system. 

Equation  (51)  is  discretized  as  follows 


[3  (E  3<»/3x)/3xl1 


*"1+1/2  *i+l/2  ~  *"1-1/2  *i-l/2 
(1/2) (x^  -  x^) 


[3  (E  34>/3y)/3y] .  .  - 

1  *  J 


*"1+1/2  *j+l/2  “  *"1-1/2  *j-l/2 
(1/2)  (yj+1  -  y^) 


[3S/3xli>j 


S1+1,J-  Sl-l,1 

xi+l  +  Xi-1 


where 


*"1+1/2  5  (L/2)  (*"i+l,1  +  ^i.j* 

*"1+1/2  2  (l/,2)  (*"i,j+l  +  £i,1^ 

*1+1/ 2  2  (*i+i,j  “  *i,j)/(Xi+l  *  Xi) 


(52) 

(53) 

(54) 


(55) 

(56) 

(57) 
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m 


*j+l/2  2  (*i,j+l  "  *i,J)/(yj+l  ~  yJJ 


(58) 


The  above  expressions  can  be  evaluated  only  for2<i<NX-l  and  for 
2  <  j  <  NY  -  1,  leaving  (NX-2)  •  (NY-2)  equations  In  NX  •  NY-4  unknowns 
(note  that  the  corner  points  of  the  grid  do  not  appear  In  the  above 
equations).  The  missing  2(NX+NY)-8  equations  are  derived  from  the  boundary 
conditions  Imposed  on  4.  For  Instance,  the  simplest  boundary  conditions 
that  could  be  Imposed  would  be  Dirichlet,  i.e.,  specification  of  known 
values  of  $  for  the  2(NX+NY)-8  grid  points  comprising  the  perimeter  of  our 
grid.  Another  possibility  would  be  Neumann  boundary  conditions,  which 
would  specify  known  values  of  the  normal  derivative  of  $  at  the  boundary. 
For  instance,  the  equation 


~  *N— l,j)/(xN  ~  *N-15  “  BXN-l/2,j  (59) 

can  be  thought  of  as  imposing  the  condition  that  the  normal  derivative  at 
the  right  boundary  point  x  ■  (1/2)  (xN  +  xN_^),  y  ■  yj  be  equal  to 
®XN-l/2  j»  the  va^ue  which  is  presumably  given. 

The  solution  to  this  linear  system  of  equations,  while  by  no  means  a 
trivial  exercise,  can  be  accomplished  by  a  number  of  algorithms.  We  have 
found  one  and  only  one  algorithm  which  will  yield  a  solution  in  all  cases, 
the  Stabilized  Error  Vector  Propagation  (SEVP)  algorithm  of  Madala 
(1978).  This  Is  a  direct  solver  and  can  be  expensive  on  a  large  grid. 
Iterative  solvers  with  which  we  have  had  success  include  the  Chebyshev 
semi-iterative  method  of  McDonald  (1980),  and  the  vectorized  incomplete 
Cholesky  conjugate  gradient  (ICCG)  algorithm  of  Hain  (1980),  which  is  an 
extension  of  the  work  of  Kershaw  (1978). 
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8.  Th«  Numerical  Solution  of  the  Continuity  Equation 

The  continuity  equation  is  ubiquitous  in  all  of  physics •  It  is  simply 
a  statement  of  the  fact  that  a  conserved  quantity  (mass  for  instance)  can 
only  appear  somewhere  in  space  if  it  comes  from  somewhere  else.  As  we  have 
noted  previously,  Eq.  (SO)  is  distinguished  by  the  appearance  of  both 
spatial  and  temporal  derivatives.  We  have  also  noted  previously  that  we 
Intend  to  treat  these  spatial  and  temporal  derivatives  in  distinctly 
different  ways  numerically.  The  formal  distinction  of  spatial  and  temporal 
derivatives  takes  the  form  of  a  general  numerical  technique  which  has  come 
to  be  known  as  the  Method  of  Lines  (MOL).  In  the  Method  of  Lines  one 
simply  treats  the  entire  spatial  differential  operator  as  some  nonlinear 
operator  H  operating  on  the  operand  or  operands  of  the  temporal  derivative 
operator,  this  case  n: 

It*  -  H(n)  (60) 


where 


H(n)  -  VA  •  (n  v^)  (61) 

Note  that  v^  is  a  function  of  n  by  Eq.  (51)  and  the  definitions 
of  Ep,  S,  and  H  is  therefore  a  very  complicated  nonlinear  operator 

acting  on  n  which  involves  all  of  the  spatial  discretization  and 
definitions  implicit  in  solving  Eq.  (51),  as  well  (as  we  shall  see)  as  the 
spatial  finite  difference  discretization  needed  to  represent  the 
operator  7^  for  Eq.  (50).  Nonetheless  this  formalism  considerably 
simplifies  our  task,  for  it  allows  us  to  cleanly  separate  out  our  treatment 
of  the  temporal  derivatives.  We  note  that  now  Eq.  (60)  is  simply  an 
ordinary  differential  equation  (ODE)  for  which  a  wide  variety  of  numerical 
Integration  techniques,  known  as  "ODE  solvers",  exist.  Actually,  as  we 
shall  see  later  Eq.  (60)  and  (61)  actually  represents  a  system  of  ODE's, 
one  for  each  spatial  grid  point,  which  are  coupled  to  each  other  through 
spatial  finite  differences  and  through  the  solution  of  the  elliptic 
equation  (51).  We  are  fortunate  here  in  that  our  system  of  ODE's  never 
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becomes  stiff  (i.e.,  there  are  no  solutions  with  time  scales  much  shorter 
than  those  of  physical  interest),  and  hence  we  have  no  need  of  the  more 
sophisticated  numerical  techniques  available  for  such  situations.  The 
solvers  actually  in  use  in  the  present  versions  of  our  simulation  codes  are 


as  follows: 

Leapfrog  -  Trapezoidal: 

n'  -  n(t-At)  +  2H(n(t)) )  At  (62a) 

n(t+At)  -  n(t)  +  (1/2)  (R(n')  +  H(n(t)))At  (62b) 

Modified  Euler: 

n'  ■  n(t)  +H(n(t))At  (63a) 

n(t+At)  -  n(t)  +  (1/2)  (H(n')  +  H(n(t)))At  (63b) 


Note  that  each  of  these  schemes  consist  of  a  predictor  (62a,  63a)  followed 
by  a  corrector  (62b,  63b),  and  that  the  corrector  stages  are  identical. 
Both  schemes  are  of  second  order  accuracy,  meaning  that  if 
n(t),  n(t-At),  H(n(t) ) ,  and  H(n(t-At))  are  known  exactly  then  the 
error  E(t+At)  in  the  solution  n(t+At)  decays  as  some  constant 
times  At2  as  At  *  0: 

E(t+At)  ♦  C  At2,  At  *  0;  C  ■  constant  (64) 

Restating  this  in  the  so-called  0  -  notation: 

E(t+At)  -  0( At2)  (65) 

The  advantage  that  the  modified  Euler  scheme  enjoys  is  that  only  n(t)  need 
be  known  to  advance  the  solution  to  time  t+At,  while  the  leapfrog- 
trapezoidal  scheme  requires  in  addition  a  knowledge  of  n(t-At).  However 
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this  advantage  la  outweighed  by  the  fact  that  the  Modified  Euler  echeae  la 
actually  weakly  unatable  for  the  caae  n(t)  -  ei0,  6  a  real  nuaber,  H(n)  • 
ikn,  k  a  poaltlve  real  nuaber.  That  la,  |n(t+At)  I  >  1  whereas 
analytically  |n(t+At)l  «1  for  all  At.  Thla  fora  of  H(n)  la  of  great 
lntereat  alnce  If  we  aet  n(x,t)  •  #i(kx-wt)  tjj#n  the  convective  derivative 
for  unit  velocity  la  3n/3x  -  lkn.  For  the  continuity  equation  thla 
lnatablllty  haa  the  effect  of  aapllfylng  the  ahort  spatial  wavelength 
coaponents  of  the  denalty  field  slightly*  The  leapfrog-trapezoidal  scheme 
does  not  have  this  defect,  and  la  therefore  the  one  we  have  chosen  for  use 
In  our  simulation  codes.  The  aodlfled  Euler  scheme  Is  used  In  our  codes 
only  to  start  the  calculation  from  the  initial  conditions,  or  to  change  the 
tiae  step,  which  must  be  done  occasionally,  alnce  even  the  leapfrog- 
trapezoidal  scheme  Is  stable  only  when  At  <  At®*n,  where  At15*11  depends  on 
the  effective  value  of  k  produced  by  the  spatial  operator  H. 

Our  problem  has  now  been  reduced  to  that  of  evaluating  the  spatial 
operator  H  on  the  finite  difference  grid  shown  In  Fig.  7.  First  we  note 
that 


H(n) 

a 

•  V±  •  (n  v^  )  -  3f/3x  +  3g/3y 

(66) 

a 

f(n) 

-  n  vx(n) 

(67) 

g(n) 

-  n  vy(n) 

(68) 

v  .  “ 
~el 

vx*+  vy  f 

(69) 

and  la  given  by  Eq.  (46).  As  we  stated  earlier,  n  and  4)  are  given  on 
the  aesh  points  (x^,  yj)  and  are  denoted  by  njj  and  i^j  respectively.  We 
shall  also  evaluate  vx  and  v^  on  these  same  grid  points,  using  centered 
finite  difference  formulae  to  be  given  presently.  Therefore  the  quantities 
f  and  g  above  are  also  known  on  these  grid  points.  We  shall  assume  for  the 
aoaant  that  our  aesh  is  uniform,  l.e.,  that  Axi+1/2  s  -  x*  is 
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independent  of  i  and  that  Ayj+1/2  =  yj+i  ~  yj  is  independent  of  j,  and 
denote  these  grid  spacings  by  simply  Ax  and  Ay  respectively*  Modifications 
necessary  for  a  nonunifora  mesh  will  be  given  later.  Then  we  can 
approximate  the  quantity  3f/3x  to  various  orders  of  accuracty: 

t||)13  •  <£i«,j  -  £u>'“  +  «“>  <70> 

*  “f11'2)  <7£> 

2££i+l,j  "  £i-l,  J7/£3fix7  ‘  ££i+2,J  “  £1-2, (72) 

+  0(Ax4) 

Similar  expressions  exist  for  approximating  3g/3y.  For  Instance 

"  (8i,j+l~  8i,j“l  )/(2Ay)  +  0(Ay2 )  (73) 

Recall  that  earlier  we  had  assumed  that  vx  and  v^  were  known  on  grid  points 
(xi,  yj).  Looking  at  Eq.  (46)  and  (47)  we  see  that  this  requires  a 

knowledge  of  7 j_$  ■  3$/3x  R  +  3 <p/3y  ?  on  grid  points,  which  are  obtained 
using  the  above  formulae  by  substituing  $  for  both  f  and  g. 

If  we  simply  choose  an  order  of  accuracy  desired  or  required  for  our 
problem,  we  have  apparently  completely  specified  our  solution  algorithm; 
and  indeed,  for  many  kinds  of  problems  this  would  be  completely 

sufficient.  However,  if  one  attempts  to  solve  even  the  simplest  of 
continuity  equations  (3n/3y  ■  0,  Vy  ■  0,  vx  ■  constant)  in  the  presence  of 

very  steep  gradients  of  n  in  the  x  direction,  the  numerical  solution  is 

soon  seen  to  be  contaminated  by  the  appearance  of  spurious  nonphysical 
oscillations  or  "ripples"  which  can  grow  in  time  and  eventually  destroy  all 
of  the  information  content  of  the  calculation.  The  reasons  are  many  and 
varied,  but  in  the  final  analysis  are  directly  caused  by  the  error 


22 


associated  with  the  finiteness  of  Ax,  Ay,  and  At:  the  "discretization 
error".  Often  this  error  can  be  reduced  to  acceptable  levels  simply  by 
using  formally  more  accurate  finite  difference  approximations  for  spatial 
and  temporal  derivatives,  for  instance  by  using  Eq.  (72)  instead  of  Eq. 
(71),  or  using  a  fourth-order  Runge  Kutta  solver  to  integrate  in  time 
instead  of  our  leapfrog-trapezoidal  scheme.  However,  if  the  spatial  or 
temporal  gradients  are  such  that  the  Taylor  series  expansion  implicit  in 
all  finite  difference  formulae  is  either  nonconvergent  or  slowly 
convergent,  then  this  technique  will  not  improve  matters  appreciably,  and 
may  even  increase  the  error.  The  brute  force  approach,  of  course,  is  to 
keep  reducing  Ax,  Ay  and  At  until  we  resolve  all  spatial  and  temporal 
structure  sufficiently  well  to  get  a  convergent  solution.  However  there 
are  many  physical  systems,  among  them  the  barium  cloud  and  equatorial 
spread  F  system,  which  allow  "shock-like”  solutions,  i.e.,  solutions  which 
contain  regions  where  the  gradient  scale  length  is  orders  of  magnitude 
smaller  that  that  of  the  other  features  in  the  problem.  On  the  scale  of 
the  overall  structure  of  the  solution,  these  regions  are  well  approximated 
by  discontinuities.  These  discontinuities  effect  the  rest  of  solution  in 
time  solely  through  their  propagation  speed  ("shock  speed”)  and  the  change 
in  the  physical  characteristics  and  velocity  of  the  plasma  across  the 
discontinuity  ("jump  conditions").  It  is  obviously  impossible  in  a 
situation  like  this  to  reduce  Ax,  Ay  and  At  to  the  point  where  the  actual 
structure  inside  the  shock  is  resolved  on  our  finite  difference  mesh. 
Fortunately,  it  is  also  unnecessary.  In  their  classic  paper,  Lax  and 
Wendroff  (1960)  showed  that  when  these  shock-like  solutions  appeared  within 
the  context  of  a  system  of  conservation  laws  (mass,  momentum,  and  energy, 
for  example),  then  any  finite  difference  scheme  which  could  represent  the 
shock  as  a  stable  propagating  entity,  regardless  of  the  computed  internal 
shock  structure,  would  recover  the  correct  shock  speed  and  jump  conditions 
(and  thus  the  correct  influence  of  the  shock  on  the  rest  of  the  solution) 
if  it  were  in  conservation  form,  a  term  we  shall  define  momentarily.  Thus 
it  is  sufficient  to  utilize  a  scheme  which  is  both  in  conservation  form  and 
which  has  the  property  of  representing  a  shock  as  a  stable  propagating 
entity.  Within  this  class  of  schemes  one  is  usually  confronted  with  a 
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choice  beCveen  schemes  which  allow  numerical  oscillations  near  the  shock 
front,  which  may  be  severe  and  which  may  In  fact  destroy  the  accuracy  of 
the  entire  calculation  if  not  carefully  controlled,  and  schemes  which 
artlfically  smear  the  shock  front  over  large  numbers  of  grid  points.  The 
oscillatory  schemes  in  general  are  of  second  or  higher  order  accuracy  In 
time  or  space,  while  the  dissipative,  non-oscillatory  schemes  are  all  first 
order  accurate  in  time  and  space.  We  shall  therefore  use  the  terms  "high 
order”  and  "low  order”  to  describe  the  above  oscillatory  and  non- 
osclllatory  schemes  respectively.  The  choice  between  high  and  low  order 
schemes  is  a  particularly  unpleasant  one.  The  inherently  high  numerical 
dissipation  of  the  low  order  schemes  tends  to  excessively  smooth  the  other 
physical  structures  in  the  problem  as  well  as  the  shock  front,  and  the  low 
convergence  rate  (0(Ax,  Ay,  At))  may  mean  that  almost  as  many  grid  points 
may  be  required  for  sufficient  accuracy  as  would  have  been  required  to 
actually  resolve  the  shock  structure  to  begin  with.  On  the  other  hand  the 
numerical  oscillations  associated  with  the  high  order  schemes  often 
propagate  into  the  entire  domain  of  the  solution,  destroying  all  of  the 
accuracy  of  the  calculation.  Again  we  are  fortunate  in  that  we  do  not  have 
to  make  this  choice,  since  we  can  have  the  best  of  both  schemes  by 
utilizing  a  technique  known  as  flux-corrected  transport  (FCT),  which  was 
originally  developed  by  Boris  and  Book  (1973)  and  later  generalized  by 
Zalesak  (1979). 

Consider  our  continuity  equation 

3n/3t  +  3f(n)/3x  +  3g(n)/3y  =■  0  (74) 

We  shall  say  that  a  finite  difference  approximation  to  Eq.  (74)  is  in 
conservation  (or  "flux”)  form  when  it  can  be  written  in  the  form 

(75) 

ni j(t+At)  "  nij(t)  ~  A  Vij  [Fi+(l/2),j'  Fi-(l/2),j+  Gi,j+(l/2) 

"  Gi » j~l/2) ] 
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where  AV^j  -  Ax^  ly  j  is  an  area  element  centered  on  grid  point  (x^, 
yj).  The  Fi+(l/2),j  an^  **i,j+(l/2)  *re  called  tranaportive  fluxes  and  are 
functions  of  f  and  g  respectively  at  one  or  aore  of  the  time  levels.  The 
functional  dependence  of  F  on  f  and  of  G  on  g  defines  the  numerical 
scheme.  For  Instance,  If  we  choose  the  trapezoidal  corrector  step  Eq. 
(62b)  combined  with  fourth  order  accurate  spatial  derivatives  Eq.  (72)  then 

Fi+(l/2),j  “  [lT  (fi+l,J  +  fljJ  “  II  (fi+2, j  +  fl-l,j)lAyj  (76) 

fi’j  *  T  (f(n"ij)  +  f(nij(t)))  (77) 

The  essence  of  the  FCT  method  is  as  follows.  For  each  time  step  one 
computes  two  distinct  sets  of  F  and  G:  one  set  by  a  low  order  scheme  (the 
"low  order  fluxes")  and  the  other  set  by  a  high  order  scheme  (the  "high 
order  fluxes").  Then  at  each  cell  Interface  (i  +  (1/2), j) 

and  (i,  j+(l/2))  one  uses  a  weighted  average  of  the  high  and  low  order  flux 
as  the  final  flux.  This  weighting  is  done  in  a  manner  which  insures  that 
the  high  order  flux  is  utilized  to  the  greatest  extent  possible  without 
introducing  numerical  oscillations  in  the  solution.  The  solution  which 
would  have  resulted  if  the  low  order  scheme  had  been  used  alone  is  used  as 
the  standard  by  which  to  determine  whether  an  oscillation  is  numerical  or 
physical.  The  result  is  a  family  of  schemes  capable  of  resolving 
discontinuities  over  2-3  grid  points  with  very  little  smearing  of  other 
physical  details  and  no  numerical  oscillations.  For  more  details,  see 

Boris  and  Book  (1973)  or  Zalesak  (1979). 

Before  closing  this  section,  let  us  briefly  describe  our  treatment  of 
nonuniform  spatial  grids.  The  basic  technique  is  to  utilize  a  smooth 
mapping  from  our  "grid  space”  (i,j)  to  real  space  (x,y).  The  mappings  we 
use  are  especially  simple  in  that  x  *  x(i)  and  y  »  y ( j ) .  Since  our 

nonuniform  spatial  mesh  enters  only  in  our  evaluation  of  3f/3x  and  g/3y  and 
since  the  treatment  for  each  is  the  same  we  shall  simply  show  our 

evaluation  of  3f/3x  here.  Utilizing  the  dummy  index  k  and  treating  it  as  a 

continuous  variable,  we  simply  use  the  chain  rule: 


Ij 


f— ) 

L3kJ 


IJ 


r _ ^ _ ) 

1  3*00^ 


rill  flif1 

l3kj1:J  l3kJt 


(78) 


Now  f  as  a  function  of  the  index  k  is  by  definition  given  on  a  uniform 
mesh,  so  we  can  use  all  of  our  previously  given  formulae  for  spatial 
derivatives  and  fluxes.  The  derivative  (3x/3k)j  can  be  taken  analytically 
if  we  have  specified  an  analytic  map  from  "k-space"  to  "x-space",  or  if 
this  map  is  not  given  explicity  but  is  still  smooth  we  can  again  use  the 
previously  given  formulae  for  spatial  derivatives  since  x  as  a  function  of 
k  is  also  by  definition  given  on  a  uniform  mesh.  In  terms  of  our  flux 
formulation,  this  simply  means  that  Ax^  is  defined  to  be  (3x/3k)^,  and  the 
rest  of  the  scheme  remains  intact. 


9.  Concluding  Remarks 

We  hope  to  have  given  the  reader  an  understanding  of  the  basic  physics 
of  the  plasma  instabilities  underlying  the  ionospheric  irregularities 
treated  here,  as  well  as  of  some  of  the  fundamental  concepts  involved  in 
the  numerical  Integration  of  the  differential  equations  describing  this 
physics.  We  cannot  treat  the  subject  in  its  entirety  here,  but  have  rather 
tried  to  give  the  reader  enough  information  to  get  started  on  his  own  if  he 
so  desires.  Both  aspects  of  the  subject,  the  physics  and  the  numerical 
analysis,  are  extremely  dynamic  fields.  Of  particular  interest  to  this 
author  is  the  fact  that  the  subject  of  numerical  solutions  to  continuity 
equations  has  recently  become  an  area  of  widespread  intensive  study  by  many 
researchers.  The  reader  is  strongly  advised  to  monitor  the  relevant 
numerical  and  mathematical  literature. 
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R-889 

Figure  1.  Photograph  of  the  Spruce  barium  cloud  24  minutes  after 
release.  Bright  areas  are  ionized  barium.  The  line  of  sight 
near  the  center  of  the  picture  is  parallel  to  the  magnetic 
field  at  the  cloud  altitude. 


27 


(ANMMMOfNOTKKp) 

t  -  0  t>0  ■■■■  * 

Figure  2.  Sketch  of  the  time  evolution  of  a  typical  barium  cloud  in  a 
plane  perpendicular  to  the  magnetic  field,  subject  to  an  upward 
directed  neutral  wind.  Lines  demarking  the  cloud  denote  plasma 
density  contours,  with  the  highest  plasma  density  in  the  center 
of  cloud. 
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Figure  3.  Two  sequential  maps  of  3  meter  radar  backscatter  at  the  earth's 
magnetic  equator.  Regions  of  Intense  backscatter  have  been 
shown  to  be  associated  with  regions  of  severe  electron  density 
depletion.  From  R.T.  Tsunoda,  J.  Geophys.  Res.,  86,  139,  1981, 
copyrighted  by  the  American  Geophysical  Union. 
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Figure  4.  Four  sequential  plots  of  contours  of  electron  density  at  the 
earth's  equator  depicting  the  formation  and  subsequent  buoyant 
rise  of  an  ESF  "bubble",  taken  from  a  numerical  simulation  of 
Zalesak  et  al.  (1982). 
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Figure  5 


.  Contours  of  constant  plasma  density  n  for  In  an  x-y  plane 
perpendicular  to  the  magnetic  field.  n  Is  a  function  of  y 
which  maximizes  at  y  ■  yQ,  superimposed  on  which  is  a 
perturbation  of  the  form  sin  kx,  where  k  is  a  wavenumber. 
Either  a  downward-directed  gravity  or  a  downward-directed 
neutral  wind  causes  ions  to  shift  slightly  leftward  relative  to 
the  electrons,  which  results  in  an  gpX  g  velocity  which  either 
damps  or  enhances  the  perturbation,  depending  on  the  sign 
of  3n/3y. 
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Figure  6.  Model  of  plasma  and  magnetic  field  geometry  used  in  this 
paper.  Field  lines  terminate  on  insulators  at  z  ■  ».  Plasma 

is  divided  into  "layers'*  along  z  for  mathematical  and  numerical 
treatment.  Each  layer  consists  of  a  single  ion  species  and  its 
associated  electrons.  Multiple  collocated  ion  species  are 
treated  by  having  multiple  collocated  "layers". 
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Spatial  finite  difference  grid  used  in  the  numerical  simulation 
codes,  showing  the  correspondence  between  i  and  x  and  between  j 
and  y.  Several  layers  of  plasma  are  shown,  even  through  the 
discussion  in  the  text  assumes  only  one  layer.  Grid  points  are 
shown  as  dots. 
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WASHINGTON,  D.C.  20360 

OICY  ATTN  NAVALEX  034  T.  HUGHES 
01CY  ATTN  PME  117 
OICY  ATTN  PME  117-T 
OICY  ATTN  CODE  5011 

COMMANDING  OFFICER 
NAVAL  INTELLIGENCE  SUPPORT  CTR 
4301  SUITLAND  ROAD,  BLDG.  5 
WASHINGTON,  D.C.  20390 

OICY  ATTN  MR.  DUBBIN  STIC  12 

OICY  ATTN  NISC-50 

OICY  ATTN  CODE  5404  J.  GALET 
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COMMANDER 

NAVAL  OCCEAN  SYSTEMS  CENTER 
SAN  DIEGO,  CA  92152 
OICY  ATTN  J.  FERGUSON 

NAVAL  RESEARCH  LABORATORY 
WASHINGTON,  D.C.  20375 

OICY  ATTN  CODE  4700  S.  L.  Ossakow 

26  CYS  IF  UNCLASS.  1  CY  IF  CLASS) 
OICY  ATTN  CODE  4701  I  Vitkovitsky 
OICY  ATTN  CODE  4780  BRANCH  HEAD  (100 
CYS  IF  UNCLASS,  1  CY  IF  CLASS) 
OICY  ATTN  CODE  7500 

OICY  ATTN  CODE  7550 

OICY  ATTN  CODE  7580 

OICY  ATTN  CODE  7551 

OICY  ATTN  CODE  7555 

OICY  ATTN  CODE  4730  E.  MCLEAN 

OICY  ATTN  CODE  4108 

OICY  ATTN  CODE  4730  B.  RIPIN 

22CY  ATTN  CODE  2628 

COMMANDER 

NAVAL  SEA  SYSTEMS  COMMAND 
WASHINGTON,  D.C.  20362 
OICY  ATTN  CAPT  R.  PITKIN 

COMMANDER 

NAVAL  SPACE  SURVEILLANCE  SYSTEM 
DAHLGREN,  VA  22448 

OICY  ATTN  CAPT  J.H.  BURTON 

OFFICER-IN-CHARGE 
NAVAL  SURFACE  WEAPONS  CENTER 
WHITE  OAK,  SILVER  SPRING,  MD  20910 
OICY  ATTN  CODE  F3i 

DIRECTOR 

STRATEGIC  SYSTEMS  PROJECT  OFFICE 
DEPARTMENT  OF  THE  NAVY 
WASHINGTON,  D.C.  20376 
OICY  ATTN  NSP-2141 
OICY  ATTN  NSSP-2722  FRED  WIMBERLY 

COMMANDER 

NAVAL  SURFACE  WEAPONS  CENTER 
DAHLGREN  LABORATORY 
DAHLGREN,  VA  22448 

OICY  ATTN  CODE  DF-14  R.  BUTLER 

OFFICER  OF  NAVAL  RESEARCH 
ARLINGTON,  VA  22217 
OICY  ATTN  CODE  465 
OICY  ATTN  CODE  461 
OICY  ATTN  CODE  402 
OICY  ATTN  CODE  420 
OICY  ATTN  CODE  421 


COMMANDER 

AEROSPACE  DEFENSE  COMMAND/DC 
DEPARTMENT  OF  THE  AIR  FORCE 
ENT  AFB,  CO  80912 

OICY  ATTN  DC  MR.  LONG 

COMMANDER 

AEROSPACE  DEFENSE  COMMAND/ XPD 
DEPARTMENT  OF  THE  AIR  FORCE 
ENT  AFB,  CO  80912 

OICY  ATTN  XPDQQ 

OICY  ATTN  XP 

AIR  FORCE  GEOPHYSICS  LABORATORY 
HANSCOM  AFB,  MA  01731 

OICY  ATTN  OPR  HAROLD  GARDNER 
OICY  ATTN  LKB  KENNETH  S.W.  CHAMPION 
OICY  ATTN  OPR  ALVA  T.  STAIR 

OICY  ATTN  PHD  JURGEN  BUCHAU 

OICY  ATTN  PHD  JOHN  P.  MULLEN 

AF  WEAPONS  LABORATORY 
KIRTLAND  AFT,  NM  87117 
OICY  ATTN  SUL 

OICY  ATTN  CA  ARTHUR  H.  GUENTHER 
OICY  ATTN  NTYCE  1LT.  G.  KRAJEI 

AFTAC 

PATRICK  AFB,  FL  32925 
OICY  ATTN  TF/MAJ  WILEY 
OICY  ATTN  TN 

AIR  FORCE  AVIONICS  LABORATORY 
WRIGHT-PATTERSON  AFB,  OH  45433 
OICY  ATTN  AAD  WADE  HUNT 
OICY  ATTN  AAD  ALLEN  JOHNSON 

DEPUTY  CHIEF  OF  STAFF 
RESEARCH,  DEVELOPMENT,  &  ACQ 
DEPARTMENT  OF  THE  AIR  FORCE 
WASHINGTON,  D.C.  20330 
OICY  ATTN  AFRDQ 

HEADQUATERS 

ELECTRONIC  SYSTEMS  DIVISION/XR 
DEPARTMENT  OF  THE  AIR  FORCE 
HANSCOM  AFB,  MA  01731 
OICY  ATTN  XR  J.  DEAS 

HEADQUATERS 

ELECTRONIC  SYSTEMS  DIVISION/YSEA 
DEPARTMENT  OF  THE  AIR  FORCE 
HANSCOM  AFB,  MA  01732 
OICY  ATTN  YSEA 
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HEADQUARTERS 

ELECTRONIC  SYSTEMS  DIVISION/DC 
DEPARTMENT  OF  THE  AIR  FORCE 
HANSCOM  AFB ,  MA  01731 

01CY  ATTN  DCKC  MAJ  J.C.  CLARK 

COMMANDER 

FOREIGN  TECHNOLOGY  DIVISION,  AFSC 
WRIGHT-PATTERSON  AFB,  OH  45433 
01CY  ATTN  NICD  LIBRARY 

01CY  ATTN  ETDP  B.  BALLARD 

COMMANDER 

ROME  AIR  DEVELOPMENT  CENTER,  AFSC 
GRIFFISS  AFB,  NY  13441 

01CY  ATTN  DOC  LIBRARY/TSLD 
01CY  ATTN  OCSE  V.  COYNE 

SAMSO/SZ 

POST  OFFICE  BOX  92960 
WORLDWAY  POSTAL  CENTER 
LOS  ANGELES,  CA  90009 
(SPACE  DEFENSE  SYSTEMS) 

01CY  ATTN  SZJ 

STRATEGIC  AIR  COMMAND/XPFS 
OFFUTT  AFB,  NB  68113 

01CY  ATTN  ADWATE  MAJ  BRUCE  BAUER 
01CY  ATTN  NRT 

01CY  ATTN  DOK  CHIEF  SCIENTIST 

SAMSO/SK 
P.O.  BOX  92960 
WORLDWAY  POSTAL  CENTER 
LOS  ANGELES,  CA  90009 

01CY  ATTN  SKA  (SPACE  COMM  SYSTEMS) 
M.  CLAVIN 

SAMSO/MN 

NORTON  AFB,  CA  92409 
(MINUTEMAN) 

01CY  ATTN  MNNL 

COMMANDER 

ROME  AIR  DEVELOPMENT  CENTER,  AFSC 
HANSCOM  AFB,  MA  01731 

OICY  ATTN  EEP  A.  LORENTZEN 

DEPARTMENT  OF  ENERGY 
LIBRARY  ROOM  G-042 
WASHINGTON,  D.C.  20545 

OICY  ATTN  DOC  CON  FOR  A.  LABOWITZ 


DEPARTMENT  OF  ENERGY 
ALBUQUERQUE  OPERATIONS  OFFICE 
P.O.  BOX  5400 
ALBUQUERQUE,  NM  87115 

OICY  ATTN  DOC  CON  FOR  D.  SHERWOOD 

EG&G,  INC. 

LOS  ALAMOS  DIVISION 

P.O.  BOX  809 

LOS  ALAMOS,  NM  85544 

OICY  ATTN  DOC  CON  FOR  J.  BREEDLOVE 

UNIVERSITY  OF  CALIFORNIA 
LAWRENCE  LIVERMORE  LABORATORY 
P.O.  BOX  808 
LIVERMORE,  CA  94550 

OICY  ATTN  DOC  CON  FOR  TECH  INFO  DEPT 
OICY  ATTN  DOC  CON  FOR  L-389  R.  OTT 
OICY  ATTN  DOC  CON  FOR  L-31  R.  HAGER 
OICY  ATTN  DOC  CON  FOR  L-46  F.  SEWARD 

LOS  ALAMOS  NATIONAL  LABORATORY 

P.O.  BOX  1663 

LOS  ALAMOS,  NM  87545 

OICY  ATTN  DOC  CON  FOR  J.  WOLCOTT 
OICY  ATTN  DOC  CON  FOR  R.F.  TASCHEK 
OICY  ATTN  DOC  CON  FOR  E.  JONES 

OICY  ATTN  DOC  CON  FOR  J.  MALIK 

OICY  ATTN  DOC  CON  FOR  R.  JEFFRIES 

OICY  ATTN  DOC  CON  FOR  J.  ZINN 

OICY  ATTN  DOC  CON  FOR  P.  KEATON 

OICY  ATTN  DOC  CON  FOR  D.  WESTERVELT 

SANDIA  LABORATORIES 
P.O.  BOX  5800 
ALBUQUERQUE,  NM  87115 

OICY  ATTN  DOC  CON  FOR  W.  BROWN 

OICY  ATTN  DOC  CON  FOR  A.  THORNBROUGH 

OICY  ATTN  DOC  CON  FOR  T.  WRIGHT 

OICY  ATTN  DOC  CON  FOR  D.  DAHLGREN 

OICY  ATTN  DOC  CON  FOR  3141 

OICY  ATTN  DOC  CON  FOR  SPACE  PROJECT  DIV 

SANDIA  LABORATORIES 
LIVERMORE  LABORATORY 
P.O.  BOX  969 
LIVERMORE,  CA  94550 

OICY  ATTN  DOC  CON  FOR  B.  MURPHEY 

OICY  ATTN  DOC  CON  FOR  T.  COOK 

OFFICE  OF  MILITARY  APPLICATION 
DEPARTMENT  OF  ENERGY 
W*.SHINGTON,  D.C.  20545 

OICY  ATTN  DOC  CON  DR.  YO  SONG 
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OTHER  GOVERNMENT  BOEING  COMPANY,  THE 

P.O.  BOX  3707 

DEPARTMENT  OF  COMMERCE  SEATTLE,  WA  98124 

NATIONAL  BUREAU  OF  STANDARDS  01CY  ATTN  G.  KEISTER 

WASHINGTON,  D.C.  20234  OICY  ATTN  D.  MURRAY 

(ALL  CORRES:  ATTN  SEC  OFFICER  FOR)  OICY  ATTN  G.  HALL 

OICY  ATTN  R.  MOORE  OICY  ATTN  J.  KENNEY 


INSTITUTE  FOR  TELECOM  SCIENCES 
NATIONAL  TELECOMMUNICATIONS  &  INFO  ADMIN 
BOULDER,  CO  80303 

OICY  ATTN  A.  JEAN  (UNCLASS  ONLY) 

OICY  ATTN  W.  UTLAUT 
OICY  ATTN  D.  CROMBIE 
OICY  ATTN  L.  BERRY 

NATIONAL  OCEANIC  &  ATMOSPHERIC  ADMIN 
ENVIRONMENTAL  RESEARCH  LABORATORIES 
DEPARTMENT  OF  COMMERCE 
BOULDER,  CO  80302 
OICY  ATTN  R.  GRUBB 
OICY  ATTN  AERONOMY  LAB  G.  REID 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

AEROSPACE  CORPORATION 
P.O.  BOX  92957 
LOS  ANGELES,  CA  90009 
OICY  ATTN  I.  GARFUNKEL 
OICY  ATTN  T.  SALMI 
OICY  ATTN  V.  JOSEPHSON 
OICY  ATTN  S.  BOWER 
OICY  ATTN  D.  OLSEN 

ANALYTICAL  SYSTEMS  ENGINEERING  CORP 
5  OLD  CONCORD  ROAD 
BURLINGTON,  MA  01803 

OICY  ATTN  RADIO  SCIENCES 

AUSTIN  RESEARCH  ASSOC.,  INC. 

1901  RUTLAND  DRIVE 
AUSTIN,  TX  78758 
OICY  ATTN  L.  SLOAN 
OICY  ATTN  R.  THOMPSON 

BERKELEY  RESEARCH  ASSOCIATES,  INC. 

P.O.  BOX  983 
BERKELEY,  CA  94701 
OICY  ATTN  J.  WORKMAN 
OICY  ATTN  C.  PRETTIE 
OICY  ATTN  S.  BRECHT 


CALIFORNIA  AT  SAN  DIEGO,  UNIV  OF 
P.O.  BOX  6049 
SAN  DIEGO,  CA  92106 
CHARLES  STARK  DRAPER  LABORATORY,  INC. 
555  TECHNOLOGY  SQUARE 
CAMBRIDGE,  MA  02139 
OICY  ATTN  D.B .  COX 
OICY  ATTN  J.P.  GILMORE 

COMSAT  LABORATORIES 
LINTHICUM  ROAD 
CLARKSBURG,  MD  20734 
OICY  ATTN  G.  HYDE 

CORNELL  UNIVERSITY 

DEPARTMENT  OF  ELECTRICAL  ENGINEERING 
ITHACA,  NY  14850 
OICY  ATTN  D.T.  FARLEY,  JR. 

ELECTROSPACE  SYSTEMS,  INC. 

BOX  1359 

RICHARDSON,  TX  75080 
OICY  ATTN  H.  LOGSTON 
OICY  ATTN  SECURITY  (PAUL  PHILLIPS) 

EOS  TECHNOLOGIES,  INC. 

606  Wlishire  Blvd. 

Santa  Monica,  Calif  90401 
OICY  ATTN  C.B.  GABBARD 

ESL,  INC. 

495  JAVA  DRIVE 
SUNNYVALE,  CA  94086 
OICY  ATTN  J.  ROBERTS 
OICY  ATTN  JAMES  MARSHALL 

GENERAL  ELECTRIC  COMPANY 
SPACE  DIVISION 
VALLEY  FORGE  SPACE  CENTER 
GODDARD  BLVD  KING  OF  PRUSSIA 
P.O.  BOX  8555 
PHILADELPHIA,  PA  19101 

OICY  ATTN  M.H.  BORTNER  SPACE  SCI  LAB 


GENERAL  ELECTRIC  COMPANY 
P.O.  BOX  1122 
SYRACUSE,  NY  13201 
OICY  ATTN  F.  REIBERT 
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GENERAL  ELECTRIC  TECH  SERVICES  CO.,  INC. 
HMES 

COURT  STREET 
SYRACUSE,  NY  13201 
01CY  ATTN  G.  HILLMAN 

GEOPHYSICAL  INSTITUTE 
UNIVERSITY  OF  ALASKA 
FAIRBANKS,  AK  99701 

(ALL  CLASS  ATTN:  SECURITY  OFFICER) 
01CY  ATTN  T.N.  DAVIS  (UNCLASS  ONLY) 
01CY  ATTN  TECHNICAL  LIBRARY 
OICY  ATTN  NEAL  BROWN  (UNCLASS  ONLY) 

GTE  SYLVANIA,  INC. 

ELECTRONICS  SYSTEMS  GRP-EASTERN  DIV 
77  A  STREET 
NEEDHAM,  MA  02194 

OICY  ATTN  DICK  STEINKOF 

HSS,  INC. 

2  ALFRED  CIRCLE 
BEDFORD,  MA  01730 

OICY  ATTN  DONALD  HANSEN 

ILLINOIS,  UNIVERSITY  OF 
107  COBLE  HALL 
150  DAVENPORT  HOUSE 
CHAMPAIGN,  IL  61820 

(ALL  CORRES  ATTN  DAN  MCCLELLAND) 

OICY  ATTN  K.  YEH 

INSTITUTE  FOR  DEFENSE  ANALYSES 
1801  NO.  BEAUREGARD  STREET 
ALEXANDRIA,  VA  22311 
OICY  ATTN  J.M.  AEIN 
OICY  ATTN  ERNEST  BAUER 
OICY  ATTN  HANS  WOLFARD 
OICY  ATTN  JOEL  BENGSTON 

INTL  TEL  &  TELEGRAPH  CORPORATION 
500  WASHINGTON  AVENUE 
NUTLEY,  NJ  07110 

OICY  ATTN  TECHNICAL  LIBRARY 

JAYCOR 

11011  TORREYANA  ROAD 
P.O.  BOX  85154 
SAN  DIEGO,  CA  92138 
OICY  ATTN  J.L.  SPERLING 


JOHNS  HOPKINS  UNIVERSITY 
APPLIED  PHYSICS  LABORATORY 
JOHNS  HOPKINS  ROAD 
LAUREL,  MD  20810 

OICY  ATTN  DOCUMENT  LIBRARIAN 
OICY  ATTN  THOMAS  POTEMRA 
OICY  ATTN  JOHN  DASSOULAS 

KAMAN  SCIENCES  CORP 
P.O.  BOX  7463 

COLORADO  SPRINGS,  CO  80933 
OICY  ATTN  T.  MEAGHER 

KAMAN  TEMPO-CENTER  FOR  ADVANCED  STUDIES 
816  STATE  STREET  (P.O  DRAWER  QQ) 

SANTA  BARBARA,  CA  93102 
OICY  ATTN  DAS I AC 
OICY  ATTN  WARREN  S.  KNAPP 
OICY  ATTN  WILLIAM  MCNAMARA 
OICY  ATTN  B.  GAMBILL 

LINKA8IT  CORP 
10453  ROSELLE 
SAN  DIEGO,  CA  92121 
OICY  ATTN  IRWIN  JACOBS 

LOCKHEED  MISSILES  &  SPACE  CO.,  INC 
P.O.  BOX  504 
SUNNYVALE,  CA  94088 
OICY  ATTN  DEPT  60-12 
OICY  ATTN  D.R.  CHURCHILL 

LOCKHEED  MISSILES  &  SPACE  CO.,  INC. 

3251  HANOVER  STREET 
PALO  ALTO,  CA  94304 

OICY  ATTN  MARTIN  WALT  DEPT  52-12 
OICY  ATTN  W.L.  IMHOF  DEPT  52-12 
OICY  ATTN  RICHARD  G.  JOHNSON  DEPT  52-12 
OICY  ATTN  J.B.  CLADIS  DEPT  52-12 

MARTIN  MARIETTA  CORP 
ORLANDO  DIVISION 
P.O.  BOX  5837 
ORLANDO,  FL  32805 
OICY  ATTN  R.  HEFFNER 


M.I.T.  LINCOLN  LABORATORY 
P.O.  BOX  73 
LEXINGTON,  MA  02173 

OICY  ATTN  DAVID  M.  TOWLE 
OICY  ATTN  L.  LOUGHLIN 
OICY  ATTN  D.  CLARK 
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MCDONNEL  DOUGLAS  CORPORATION 
5301  BOLSA  AVENUE 
HUNTINGTON  BEACH,  CA  92647 
OICY  ATTN  N.  HARRIS 
OICY  ATTN  J.  HOULE 
01CY  ATTN  GEORGE  MtOZ 
OICY  ATTN  H.  OLSON 
OICY  ATTN  R.W.  HAL PR IN 
OICY  ATTN  TECHNICAL  LIBRARY  SERVICES 

MISSION  RESEARCH  CORPORATION 
735  STATE  STREET 
SANTA  BARBARA,  CA  93101 
OICY  ATTN  P.  FISCHER 
OICY  ATTN  W.F.  CREVIER 
OICY  ATTN  STEVEN  L.  GUTSCHE 
OICY  ATTN  D.  SAPPENFIELD 
OICY  ATTN  R.  BOGUSCH 
OICY  ATTN  R.  HENDRICK 
OICY  ATTN  RALPH  KIL8 
OICY  ATTN  DAVE  SOWLE 
OICY  ATTN  F.  FAJEN 
OICY  ATTN  M.  SCHEIBE 
OICY  ATTN  CONRAD  L.  LONGMIRE 
OICY  ATTN  B.  WHITE 

MISSION  RESEARCH  CORP. 

1400  SAN  MATEO  BLVD.  SE 
SUITE  A 

ALBUQUERQUE,  NEW  MEXICO  87108 
OICY  R.  STELLINGWERF 
OICY  M.  ALME 
OICY  L.  WRIGHT 

MITRE  CORPORATION,  THE 
P.O.  BOX  208 
BEDFORD,  MA  01730 

OICY  ATTN  JOHN  MORGANSTERN 
OICY  ATTN  G.  HARDING 
OICY  ATTN  C.E.  CALLAHAN 

MITRE  CORP 

WESTGATE  RESEARCH  PARK 
1820  DOLLY  MADISON  BLVD 
MCLEAN,  VA  22101 
OICY  ATTN  W.  HALL 
OICY  ATTN  W.  FOSTER 

PACIFIC-SIERRA  RESEARCH  CORP 
12340  SANTA  MONICA  BLVD. 

LOS  ANGELES,  CA  90025 

OICY  ATTN  E.C.  FIELD,  JR. 


PENNSYLVANIA  STATE  UNIVERSITY 
IONOSPHERE  RESEARCH  LAB 
318  ELECTRICAL  ENGINEERING  EAST 
UNIVERSITY  PARK,  PA  16802 
(NO  CLASS  TO  THIS  ADDRESS) 

OICY  ATTN  IONOSPHERIC  RESEARCH  LAB 

PHOTOMETRICS ,  INC. 

4  ARROW  DRIVE 
WOBURN,  MA  01801 

OICY  ATTN  IRVING  L.  KOFSKY 

PHYSICAL  DYNAMICS,  INC. 

P.O.  BOX  3027 
BELLEVUE,  WA  98009 

OICY  ATTN  E.J.  FREMOUW 

PHYSICAL  DYNAMICS,  INC. 

P.O.  BOX  10367 
OAKLAND,  CA  94610 
ATTN  A.  THOMSON 

R  &  D  ASSOCIATES 
P.O.  BOX  9695 
MARINA  DEL  RET,  CA  90291 
OICY  ATTN  FORREST  GILMORE 
OICY  ATTN  WILLIAM  B.  WRIGHT,  JR. 
OICY  ATTN  ROBERT  F.  LELEVIER 
OICY  ATTN  WILLIAM  J.  KARZAS 
OICY  ATTN  H.  ORY 
OICY  ATTN  C.  MACDONALD 
OICY  ATTN  R.  TURCO 
OICY  ATTN  L.  DcRAND 
OICY  ATTN  W.  TSAI 

RAND  CORPORATION,  THE 
1700  MAIN  STREET 
SANTA  MONICA,  CA  90406 
OICY  ATTN  CULLEN  CRAIN 
OICY  ATTN  ED  BEDROZIAN 

RAYTHEON  CO. 

528  BOSTON  POST  ROAD 
SUDBURY,  MA  01776 

OICY  ATTN  BARBARA  ADAMS 

RIVERSIDE  RESEARCH  INSTITUTE 
80  WEST  END  AVENUE 
NEW  YORK,  NY  10023 

OICY  ATTN  VINCE  TRAPANI 
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SCIENCE  APPLICATIONS,  INC. 
P.O.  BOX  23SI 
LA  JOLLA,  CA  92038 
01CT  ATTN  LEWIS  M.  LINSON 
01CT  ATTN  DANIEL  A.  HAMLIN 
Old  ATTN  E.  PRIEMAN 
Old  ATTN  E.A.  STRAKER 
Old  ATTN  CURTIS  A.  SMITH 
Old  ATTN  JACK  MCDOUGALL 

SCIENCE  APPLICATIONS,  INC 
1710  GOODRIDGE  DR. 

MCLEAN,  VA  22102 
ATTN:  J.  COCKAYNE 

SRI  INTERNATIONAL 
333  RAVENSNOOD  AVENUE 
MENLO  PARK,  CA  94025 


Old 

ATTN 

DONALD  NEILSON 

Old 

ATTN 

ALAN  BURNS 

Old 

ATTN 

G.  SMITH 

Old 

ATTN 

R.  TSUNODA 

Old 

ATTN 

DAVID  A.  JOHNSON 

Old 

ATTN 

WALTER  G.  CHESNUT 

Old 

ATTN 

CHARLES  L.  RINO 

Old 

ATTN 

WALTER  JAYE 

Old 

ATTN 

J.  VICKREY 

Old 

ATTN 

RAY  L.  LEAD AB RAND 

Old 

ATTN 

G.  CARPENTER 

Old 

ATTN 

G.  PRICE 

Old 

ATTN 

J.  PETERSON 

Old 

ATTN 

R.  LIVINGSTON 

01CY 

ATTN 

V.  GONZALES 

Old 

ATTN 

D.  MCDANIEL 

STEWART  RADIANCE  LABORATORY 
UTAH  STATE  UNIVERSITY 
1  DE  ANGELO  DRIVE 
BEDFORD,  MA  01730 
Old  ATTN  J.  ULWICK 

TECHNOLOGY  INTERNATIONAL  CORP 
75  WIGGINS  AVENUE 
BEDFORD,  MA  01730 

Old  ATTN  W.P.  BOqUIST 

TOYON 

34  WALNUT  LAND 
SANTA  BARBARA,  CA  93111 
Old  ATTN  JOHN  ISE,  JR. 
Old  ATTN  JOEL  GARBARINO 


TRW  DEFENSE  A  SPACE  SYS  GROUP 
ORE  SPACE  PARK 
REDONDO  BEACH,  CA  9C’78 
Old  ATTN  R.  K.  PLEBUCH 
Old  ATTN  S.  ALTSCHULER 
Old  ATTN  0.  DEB 
Old  ATTN  D.  STOCKWELL 
SNTF/1575 

VIS1DYNE 

SOOTH  BEDFORD  STREET 
BURLINGTON,  MASS  01803 
Old  ATTN  W.  REIDY 
Old  ATTN  J.  CARPENTER 
Old  ATTN  C.  HUMPHREY 
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IONOSPHERIC  MODELING  DISTRIBUTION  LIST 
(UNCLASSIFIED  ONLY) 


PLEASE  DISTRIBUTE  ONE  COPY  TO  EACH  OF  THE 
NOTED) 

NAVAL  RESEARCH  LABORATORY 
WASHINGTON,  D.C.  20375 
DR.  P.  MANGE  -  CODE  4101 
DR.  E.  SZUSZCZEUICZ  -  CODE  4108 
DR.  J.  GOODMAN  -  CODE  4180 
DR.  P.  RODRIGUEZ  -  CODE  4108 

A.F.  GEOPHYSICS  LABORATORY 
L.G.  HANSCOM  FIELD 
BEDFORD,  MA  01730 
DR.  T.  ELKINS 
DR.  W.  SUIDER 
MRS.  R.  SAGALYN 
DR.  J.M.  FORBES 
DR.  T.J.  KENESHEA 
DR.  W.  BURKE 
DR.  H.  CARLSON 
DR.  J.  JASPERSE 

BOSTON  UNIVERSITY 
DEPARTMENT  OF  ASTRONOMY 
BOSTON,  MA  02215 
DR.  J.  AARONS 

CORNELL  UNIVERSITY 
ITHACA,  NY  14850 
DR.  W.E.  SWARTZ 
DR.  R.  SUDAN 
DR.  D.  FARLEY 
DR.  M.  KELLEY 

HARVARD  UNIVERSITY 
HARVARD  SQUARE 
CAMBRIDGE,  MA  02138 
DR.  M.B.  McELROY 
DR.  R.  LINDZEN 

INSTITUTE  FOR  DEFENSE  ANALYSIS 
400  ARMY/NAVY  DRIVE 
ARLINGTON,  VA  22202 
DR.  E.  BAUER 


FOLLOWING  PEOPLE  (UNLESS  OTHERWISE 


MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
PLASMA  FUSION  CENTER 
LIBRARY,  NW1 6-262 
CAMBRIDGE,  MA  02139 

NASA 

GODDARD  SPACE  FLIGHT  CENTER 
GREENBELT,  MD  20771 
DR.  R.F.  BENSON 
DR.  K.  MAEDA 
Dr.  S.  CURTIS 
Dr.  M.  DUB IN 

DR.  N.  MAYNARD  -  CODE  696 

NATIONAL  TECHNICAL  INFORMATION  CENTER 
CAMERON  STATION 
ALEXANDRIA,  VA  22314 
12CY  ATTN  TC 

COMMANDER 

NAVAL  AIR  SYSTEMS  COMMAND 
DEPARTMENT  OF  THE  NAVY 
WASHINGTON,  D.C.  20360 
DR.  T.  CZUBA 

COMMANDER 

NAVAL  OCEAN  SYSTEMS  CENTER 
SAN  DIEGO.  CA  92152 

MR.  R.  ROSE  -  CODE  5321 

NOAA 

DIRECTOR  OF  SPACE  AND  ENVIRONMENTAL 
LABORATORY 
BOULDER,  CO  80302 
DR.  A.  GLENN  JEAN 
DR.  G.W.  ADAMS 
DR.  D.N.  ANDERSON 
DR.  K.  DAVIES 
DR.  R.  F.  DONNELLY 

OFFICE  OF  NAVAL  RESEARCH 
800  NORTH  QUINCY  STREET 
ARLINGTON,  VA  22217 
DR.  G.  JOINER 


PENNSYLVANIA  STATE  UNIVERSITY 
UNIVERSITY  PARK,  PA  16802 
DR.  J.S.  NISBET 
DR.  P.R.  ROHRBAUGH 
DR.  L.A.  CARPENTER 
DR.  M.  LEE 
DR.  R.  DIVANY 
DR.  P.  BENNETT 
DR.  F.  KLEVANS 

PRINCETON  UNIVERSITY 
PLASMA  PHYSICS  LABORATORY 
PRINCETON,  NJ  08540 
DR.  F.  PERKINS 

SCIENCE  APPLICATIONS,  INC. 

1150  PROSPECT  PLAZA 
LA  JOLLA,  CA  92037 
DR.  D.A.  HAMLIN 
DR.  L.  LINSON 
DR.  E.  FRIEMAN 

STANFORD  UNIVERSITY 
STANFORD,  CA  94305 
DR.  P.M.  BANKS 

U.S.  ARMY  ABERDEEN  RESEARCH 
AND  DEVELOPMENT  CENTER 
BALLISTIC  RESEARCH  LABORATORY 
ABERDEEN,  MD 

DR.  J.  HEIMERL 

GEOPHYSICAL  INSTITUTE 
UNIVERSITY  OF  ALASKA 
FAIRBANKS,  AK  99701 
DR.  L.C.  LEE 

UNIVERSITY  OF  CALIFORNIA, 
BERKELEY 

BERKELEY,  CA  94720 
DR.  M.  HUDSON 

UNIVERSITY  OF  CALIFORNIA 
LOS  ALAMOS  SCIENTIFIC  LABORATORY 
J-10,  MS-664 
LOS  ALAMOS,  NM  87545 
DR.  M.  PONCRATZ 
DR.  D.  SIMONS 
DR.  G.  BARAS CH 
DR.  L.  DUNCAN 
DR.  P.  BERNHARDT 
DR.  S.P.  CARY 


UNIVERSITY  OF  CALIFORNIA, 
LOS  ANGELES 
405  HILLGARD  AVENUE 
LOS  ANGELES,  CA  90024 
DR.  F.V.  CORONITI 
DR.  C.  KENNEL 
DR.  A.Y.  WONG 

UNIVERSITY  OF  MARYLAND 
COLLEGE  PARK,  MD  20740 
DR.  K.  PAPADOPOULOS 
DR.  E.  OTT 

JOHNS  HOPKINS  UNIVERSITY 
APPLIED  PHYSICS  LABORATORY 
JOHNS  HOPKINS  ROAD 
LAUREL,  MD  20810 
DR.  R.  GREENWALD 
DR.  C.  MENG 

UNIVERSITY  OF  PITTSBURGH 
PITTSBURGH,  PA  15213 
DR.  N.  ZABUSKT 
DR.  M.  BIONDI 
DR.  E.  OVERMAN 

UNIVERSITY  OF  TEXAS 
AT  DALLAS 

CENTER  FOR  SPACE  SCIENCES 
P.O.  BOX  688 

RICHARDSON,  TEXAS  75080 
DR.  R.  HEELIS 
DR.  W.  HANSON 
dr.  j.p.  McClure 

UTAH  STATE  UNIVERSITY 
4TH  AND  8TH  STREETS 
LOGAN,  UTAH  84322 
DR.  R.  HARRIS 
DR.  K.  BAKER 
DR.  R.  S CHUNK 
DR.  J.  ST. -MAURICE 


KIRUKA  GEOPHYSICAL  INSTITUTE 
BOX  709 

S- 98 127  KIHUNA,  SWEDES 
CHRISTEB  JUREN 
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